#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#define ll long long
#define maxn 100+5
#pragma GCC optimize(2)
const double eps=1e-11;

void printMartix(double c[],double a[],double d[],int N);
void LU(double c[],double a[],double d[],double p[],double q[],double r[],int N);
void solve(double x[],double y[],double q[],double p[],double r[],double b[],int N);

int main(){
    int n;
    double a[maxn],d[maxn],c[maxn],b[maxn];
    double q[maxn],p[maxn],r[maxn];
    double x[maxn],y[maxn];
    double Zero[maxn],One[maxn];
    printf("请输入矩阵阶数n:");
    scanf("%d",&n);
    for(int i=1;i<=n;++i) Zero[i]=0.0;
    for(int i=1;i<=n;++i) One[i]=1.0;
    printf("请输入%d维向量a(主对角线)\n",n);
    for(int i=1;i<=n;++i)
        scanf("%lf",&a[i]);
    printf("请输入%d维向量d(主对角线下面一行)\n",n-1);
    for(int i=2;i<=n;++i)
        scanf("%lf",&d[i]);
    printf("请输入%d维向量c(主对角线上面一行)\n",n-1);
    for(int i=1;i<=n-1;++i)
        scanf("%lf",&c[i]);
    printf("生成的三对角矩阵为:\n");
    printMartix(c,a,d,n);
    LU(c,a,d,p,q,r,n);
    printf("Crout分解完成,其中L=\n");
    printMartix(Zero,p,r,n);
    printf("U=\n");
    printMartix(q,One,Zero,n);
    printf("输入%d维的方程右侧向量b\n",n);
    for(int i=1;i<=n;++i)
        scanf("%lf",&b[i]);
    solve(x,y,q,p,r,b,n);
    printf("方程组的解为：\n");
    for(int i=1;i<=n;++i)
        printf("x%d=%10.6lf\n",i,x[i]);
    return 0;
}

inline void printMartix(double c[],double a[],double d[],int N){
    for(int i=1;i<=N;++i){
        for(int j=1;j<=N;++j){
            if(i==j)
                printf("%10.6lf",a[i]);
            else if(i==j-1)
                printf("%10.6lf",c[i]);
            else if(i==j+1)
                printf("%10.6lf",d[i]);
            else   
                printf("%10.6lf",0.00);
        }
        printf("\n");
    }
}

inline void LU(double c[],double a[],double d[],double p[],double q[],double r[],int N){
    p[1]=a[1];
    for(int i=2;i<=N;++i)
        r[i]=d[i];
    for(int i=1;i<=N-1;++i){
        q[i]=c[i]/p[i];
        p[i+1]=a[i+1]-d[i+1]*q[i];
    }
}

inline void solve(double x[],double y[],double q[],double p[],double r[],double b[],int N){
    y[1]=b[1]/p[1];
    for(int i=2;i<=N;++i)
        y[i]=(b[i]-r[i]*y[i-1])/p[i];
    
    x[N]=y[N];
    for(int i=N-1;i>=1;--i)
        x[i]=y[i]-q[i]*x[i+1];
}
